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ABSTRACT 

We  present  a  detailed  description  of  the  development,  implementation,  and 
application  of  a  computer  program  to  calculate  the  detonation  parameters  of 
condensed  phase  explosives.  The  code  is  based  on  Mader's  BKW  chemical  equilibrium 
code,  but  contains  important  new  features.  A  new  algorithm  to  calculate  the  minimum 
in  the  free  energy  of  the  product  composition  has  been  included.  This  is  a  probabilistic 
algorithm,  based  on  the  method  of  Benke  and  Skinner,  and  its  inclusion  ensures  that 
the  true  global  minimum  in  the  free  energy  will  always  be  found.  As  well  as  the  BKW 
equation  of  state  to  describe  the  detonation  products,  the  new  code  also  includes  the 
JCZ3  equation  of  state.  This  is  an  intermolecular  equation  of  state  containing  no 
adjustable  parameters,  and  hence  should  be  applicable  to  a  wider  range  of  explosives 
than  could  be  considered  using  the  BKW  code.  We  have  validated  the  code  on  a  wide 
range  of  military  explosives,  using  both  the  new  probabilistic  minimisation  algorithm 
as  well  as  the  original  method  of  steepest  descent,  for  both  the  BKW  and  JCZ3 
equations  of  state.  We  also  present  a  detailed  description  of  the  application  of  the  code 
to  the  non-ideal  underwater  explosive  PBXN-111,  and  show  that  the  performance  of 
the  explosive  is  best  described  using  the  JCZ3  equation  of  state. 
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Chemical  Equilibrium  Calculations  for 
Detonation  Products 


Executive  Summary 


This  report  describes  the  development  and  application  of  a  chemical  equilibrium 
computer  code  to  calculate  the  detonation  parameters  of  condensed  phase  explosives. 
The  code  is  based  on  the  Chapman-Jouguet  theory  of  detonation,  incorporates  two 
different  equations  of  state  to  describe  the  detonation  products,  and  offers  a  choice  of 
two  different  algorithms  to  minimize  the  free  energy  of  the  product  composition. 

We  first  describe  the  development  of  the  code  using  the  Becker-Kistiakowsky-Wilson 
(BKW)  equation  of  state  to  describe  the  detonation  products,  and  the  method  of 
steepest  descent  to  minimize  the  free  energy  of  the  product  composition.  The  code  is 
then  used  to  calculate  detonation  parameters  for  a  number  of  standard  ideal  military 
explosives  and  shown  to  give  results  which  are  identical  to  those  obtained  from 
existing  codes  employing  the  same  equation  of  state  and  minimization  algorithm. 

One  weakness  of  the  method  of  steepest  descent  is  that  it  finds  local  minima,  and  to 
obtain  accurate  predictions  of  the  detonation  parameters  it  must  be  given  an  initial 
estimate  of  the  detonation  state  which  is  relatively  close  to  the  final  equilibrium  state. 
For  ideal  military  explosives  this  causes  no  problems,  but  one  of  the  objectives  in 
developing  the  current  code  was  to  apply  it  to  highly  non-ideal  underwater  explosives 
such  as  PBXN-111,  in  which  case  the  procedure  is  far  from  straightforward.  Hence 
there  is  a  finite  probability  that  the  algorithm  will  find  a  local  minimum,  rather  than 
the  global  minimum,  and  thereby  give  a  misleading  result  for  the  detonation  velocity 
and  pressure. 

To  eliminate  this  possibility  we  also  added  a  second  algorithm  to  the  code  to  minimize 
the  free  energy.  This  is  a  stochastic  method  which  uses  an  adaptive  probabilistic 
technique  to  locate  minima,  and  is  guaranteed  to  find  the  global  minimum  without  the 
need  to  make  an  initial  estimate  of  the  final  equilibrium  state.  We  checked  that  the 
new  algorithm  gave  identical  results  to  the  steepest  descent  technique  for  the  ideal 
explosives  already  considered,  and  then  used  the  code  to  calculate  the  detonation 
parameters  for  PBXN-111.  Both  algorithms  gave  essentially  the  same  result,  verifying 
that  the  true  equilibrium  composition  was  being  found,  but  the  calculated  detonation 
pressure  and  velocity  were  significantly  different  from  the  experimental  values. 

To  improve  agreement  with  experiment  we  then  implemented  the  JCZ3  equation  of 
state  in  the  code.  JCZ3  is  based  on  an  intermolecular  potential  and  contains  no 
adjustable  parameters,  and  hence  was  considered  to  be  more  appropriate  for  highly 
non-ideal  explosives  than  the  BKW  equation  of  state.  Implementation  of  JCZ3  was 
validated  on  a  variety  of  standard  explosives,  and  then  applied  to  calculate  the 
detonation  parameters  of  PBXN-111.  This  resulted  in  much  better  agreement  with 
experiment,  and  also  in  very  good  agreement  with  calculations  using  related 
equations  of  state  based  on  intermolecular  potentials. 
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1.  Introduction 

Calculations  based  on  the  Chapman-Jouguet  (CJ)  theory  of  steady  state  detonation 
have  proven  very  successful  in  the  past  in  reproducing  the  experimentally  determined 
detonation  velocities  and  pressures  of  ideal  explosive  compositions  [1,2].  The  CJ 
theory  assumes  that  the  flow  is  one-dimensional,  the  reaction  zone  is  infinitely  thin, 
and  that  completion  of  reaction  coincides  with  the  sonic  point  in  the  flow  (which  is 
then  defined  as  the  CJ  point).  Calculation  of  the  detonation  pressure  and  velocity  at 
the  CJ  point  then  requires  the  use  of  a  chemical  equilibrium  code  to  determine  the 
product  composition  with  the  minimum  free  energy  (which  we  refer  to  as  minimizing 
the  free  energy),  and  appropriate  equations  of  state  to  describe  the  detonation 
products  at  the  very  high  temperatures  and  pressures  involved  in  the  detonation 
process.  There  are  a  number  of  computer  programs  capable  of  performing  these 
calculations  and  for  many  years  AMRL  has  been  using  a  version  of  Mader's  BKW  code 
for  this  purpose  [3].  When  applied  to  ideal  CHNO  explosives  the  agreement  between 
calculated  and  experimental  detonation  velocity  and  pressure  has  been  excellent  [4]. 

Recently  in  Australia  there  has  been  considerable  interest  in  the  composite  explosive 
PBXN-111  (previously  known  as  PBXW-115)  for  use  in  underwater  munitions  [5-9]. 
This  is  a  highly  non-ideal  composition  however,  and  the  use  of  a  chemical  equilibrium 
code  to  calculate  a  detonation  velocity  and  pressure  for  such  explosives  is  far  from 
straightforward.  Mader  has  made  some  progress  in  this  area  by  using  the  BKW 
chemical  equilibrium  code  and  limiting  the  extent  of  reaction  of  some  of  the 
constituents.  This  mimics  the  effect  of  the  late-time  reactions  which  occur  in  non-ideal 
explosives,  and  Mader  has  been  able  to  get  excellent  agreement  with  experimental 
data  for  detonation  velocity  and  pressure  using  this  method  [1].  This  procedure 
determines  the  equation  of  state  of  the  detonation  products,  which  can  then  be  used  in 
a  hydrocode  calculation  to  simulate  aquarium  test  data  and  predict  the  position  of 
measured  shock  wave  and  confinement/water  interface  positions.  Good  agreement 
with  some  of  the  experimental  data  has  been  obtained  using  this  approach  [10],  but 
the  method  has  the  disadvantage  of  requiring  a  new  determination  of  the  equation  of 
state  when  the  experimental  conditions  are  changed. 

Jones  and  Kennedy  [11]  have  adopted  a  different  approach  to  the  modeling  of  highly 
non-ideal  composite  explosives  by  using  the  slightly  divergent  reactive  flow  theory  of 
Kirby  and  Leiper  [12],  together  with  the  reactive  hydrocode  DYNA2D.  This  method 
has  the  advantage  of  allowing  accurate  resolution  of  the  reaction  zone  length  in 
composite  explosives  by  assuming  a  more  realistic  expression  for  the  rate  of  energy 
release  based  on  a  sequence  of  decomposition  reactions  occurring  on  different  time 
scales.  Appropriate  values  for  the  time  constants  appearing  in  these  equations  are 
then  obtained  from  experimental  data  on  the  variation  of  detonation  velocity  with 
charge  diameter.  The  method  also  requires  an  estimate  of  the  detonation  velocity  at 
infinite  diameter,  and  this  is  obtained  from  a  chemical  equilibrium  code,  but  in  this 
application  the  need  to  restrict  the  degree  of  reaction  of  some  of  the  constituents  of  the 
equilibrium  composition  is  removed.  In  applying  this  approach  however  it  has  been 
found  that  the  BKW  code  appears  to  overestimate  the  detonation  velocity  of  some 
composite  explosives.  For  ANFO-like  explosives  for  example  Kennedy  has  found  that 
the  JCZ3  equation  of  state  in  ICI's  detonation  code  HEDEQ  gives  better  estimates  for 
the  detonation  velocity  of  these  explosives  [13],  and  Kennedy  and  Jones  [14]  have 
found  that  the  detonation  velocity  of  PBXN-111  is  more  accurately  calculated  using  the 
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ICI  chemical  equilibrium  code  IDEX,  which  uses  an  intermolecular  equation  of  state 
for  the  gaseous  products. 

Mader's  BKW  code  employs  the  Becker-Kistiakowsky-Wilson  equation  of  state  for  the 
gaseous  detonation  products  [15],  and  the  Cowan-Fickett  equation  of  state  for  the  solid 
products  [16].  The  BKW  equation  of  state  is  a  phenomenological  expression  containing 
several  adjustable  parameters  which  have  been  determined  by  fitting  to  a  series  of 
experimental  results  obtained  on  ideal  explosives.  More  realistic  equations  of  state 
based  on  intermolecular  potentials  and  refined  mixing  rules  are  now  available,  but  the 
effect  of  these  different  equations  of  state  on  the  calculated  detonation  properties  of 
highly  non-ideal  explosives  such  as  PBXN-111  has  not  been  extensively  investigated. 
The  BKW  code  is  not  easily  modified  however,  and  could  not  be  used  to  study  the 
effect  of  different  equations  of  state  on  the  detonation  performance  of  PBXN-111. 

Determination  of  the  CJ  state  of  an  explosive  requires  the  minimization  of  the  free 
energy  of  the  product  composition.  The  BKW  code  uses  the  method  of  steepest 
descent  to  perform  this  calculation  [17],  and  this  method  requires  an  initial  estimate  of 
the  CJ  state  to  be  made.  The  method  of  steepest  descent  then  finds  the  closest 
minimum  to  the  initial  guess,  but  is  not  guaranteed  to  find  the  global  minimum.  In 
actual  practice,  when  dealing  with  CHNO  explosives,  the  method  invariably  does  find 
the  global  minimum  because  realistic  initial  estimates  of  the  CJ  state  can  be  made.  For 
highly  non-ideal  explosives  this  is  not  necessarily  the  case  however,  and  a  method  of 
minimizing  the  free  energy  which  was  guaranteed  to  find  the  global  minimum 
without  the  need  to  make  an  initial  estimate  of  the  CJ  state  would  be  highly  desirable. 

The  above  considerations  have  lead  us  to  develop  our  own  chemical  equilibrium  code 
for  the  determination  of  the  CJ  state  of  condensed  phase  explosives.  The  code  is 
written  in  modular  form  and  the  equation  of  state  has  been  incorporated  as  a 
subroutine  so  that  it  can  easily  be  replaced  by  a  new  equation  of  state  by  rewriting  this 
subroutine.  Two  numerical  techniques  have  been  used  to  minimize  the  free  energy; 
we  use  the  method  of  steepest  descent  as  used  by  Mader,  both  because  of  its  speed  of 
execution,  and  the  need  to  check  agreement  with  Mader's  calculations,  and  we  also 
use  a  new  probabilistic  algorithm  of  Benke  and  Skinner  [18],  which  is  designed  to  find 
the  global  minimum  without  the  need  for  an  initial  estimate  of  the  equilibrium 
composition. 

In  Section  2  we  outline  the  equations  describing  the  thermodynamics  of  the 
equilibrium  composition  and  their  solution  using  the  method  of  steepest  descent  and 
the  BKW  equation  of  state.  In  Section  3  we  describe  the  probabilistic  algorithm  of 
Benke  and  Skinner  and  it's  implementation  into  our  chemical  equilibrium  code. 
Calculations  of  the  CJ  state  of  standard  CHNO  explosives  using  our  code  are  in 
excellent  agreement  with  Mader's  BKW  results  when  either  the  method  of  steepest 
descent  or  the  probabilistic  algorithm  are  used  for  the  free  energy  minimization.  In 
Section  4  we  use  our  code  to  calculate  the  infinite  diameter  detonation  velocity  and 
pressure  for  PBXN-111,  and  also  consider  the  effect  of  limiting  the  extent  of  reaction  of 
some  of  the  products  on  the  calculated  detonation  properties.  In  Section  5  we  describe 
the  implementation  of  the  JCZ3  intermolecular  equation  of  state  [19]  into  our  chemical 
equilibrium  code,  and  then  discuss  the  effect  which  this  has  on  the  calculated 
detonation  parameters  of  PBXN-111. 
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2.  Determining  the  CJ  State  of  Explosives 


2.1  CJ  Condition  and  Conservation  Relation 

The  CJ  state  of  an  explosive  can  be  determined  by  satisfying  the  following  criteria 
[1/2]: 

1.  At  the  CJ  point  the  slope  of  the  Raleigh  line  is  tangent  to  the  isentrope,  ie 


dP  _  P0  -P 
9V  s  “  V0  -  V 


(1) 


2.  Energy  must  be  conserved  across  the  shock  front,  ie  the  Hugoniot  condition 
must  be  satisfied. 


E-E0  =|(P  +  P„)(V0-V)  (2) 

In  order  to  apply  these  criteria  we  need  to  be  able  to  compute  the  internal  energy  of 
the  equilibrium  composition  of  the  detonation  products.  The  equilibrium  composition 
can  be  calculated  by  minimizing  the  free  energy  of  the  mixture  of  detonation  products. 


2.2  Computation  of  Chemical  Equilibrium 

For  a  complex  mixture  of  reaction  products  the  equilibrium  composition  is  determined 
by  minimizing  the  free  energy  of  the  products  within  the  constraint  of  mass  balance. 
The  free  energy  of  the  mixture  can  be  expressed  as  the  sum  of  the  free  energy  of  the 
components: 


F(X)  =  Ifj  (3) 

i=l 

where  X  represents  the  product  composition  vector.  The  free  energy  fj  is  given  by 

fj  =  Fi  +  /n[xQ  (4) 

f  * 

where  Xj  is  the  mole  fraction  of  the  ith  component.  The  expression  for  Fj  will 

depend  on  the  particular  equation  of  state  chosen  to  describe  the  detonation  products. 

To  obtain  the  equilibrium  composition,  X  =  (xj,  X2 . xN)  where  Xj  is  the  number  of 

moles  of  the  i™  product,  then  F(X)  is  minimized  subject  to  the  mass  balance 
constraint.  A  commonly  used  method  to  achieve  this  is  the  method  of  steepest 
descent,  as  described  by  White  et  al.  [17].  We  start  with  any  positive  set  of  values  Y 
=(yi,  y2....yN>  which  satisfy  the  mass  balance  equation.  The  free  energy  is  then 
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expanded  in  a  Taylor's  series  to  second  order  about  Y  and  gives  Q(X),  the  quadratic 
approximation  to  F(X)/  as  follows 


Q(X)  =  F(Y)  +  Ai 


3F 


T3xi 


+  - AiAkH 


a2F 


X=Y 


k  3xidxk 


(5) 


where  Aj  =  xj  -  yj  .  To  minimize  Q(X)  subject  to  the  mass  balance  constraint,  ie. 


N 

L«ikxi=bk  <6> 

i=l 


we  form  the  function  G(X)  defined  by 


G(X)  =  Q(X)  +  Xtt 
j 


f 

— X  aijxi  “l"  b  j 
V  i 


(7) 


where  the  rcj  are  Lagrange  multipliers.  To  minimize  G(X),  we  set  d  G/d  Xj  =  0  and 
find  the  following  set  of  equations: 

— -^+  I^k«ik  =  -fi(V)  (8) 

yt  y  k=i 


N  N 

where  x=Xxj,  y=£yi  (9) 

i=l  i=l 


and  M  is  the  number  of  elements,  N  is  the  number  of  products,  Ojk  is  the  product 
elemental  composition  matrix,  and  bk  is  the  reactant  elemental  composition  vector. 
Equations  (8)  and  (6)  are  then  used  to  construct  a  matrix  equation  which  can  be  solved 
by  one  of  the  standard  methods.  We  use  the  subroutine  GAUSS  as  described  in  Press 
et  al.  [20].  The  calculation  of  the  equilibrium  composition  for  a  given  P  and  T  is  an 
iterative  procedure.  An  initial  estimate  of  the  composition  vector  (Y  =  yj,  y2--yN)  is 
made  and  then  equations  (6)  and  (8)  solved  to  provide  the  new  estimate  X  =  (xj, 
X2-..X]\j).  The  free  energy  of  the  new  composition  is  then  calculated  and  equations  (6) 
and  (8)  again  solved  with  this  new  value  to  give  a  new  estimate  of  the  composition. 
This  procedure  is  repeated  until  the  difference  between  X  and  Y  falls  below  some 
predetermined  level. 
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2.3  Thermodynamic  Functions  for  the  BKW  Equation  of  State 

In  order  to  determine  the  equilibrium  composition,  expressions  for  the  free  energy  of 
the  detonation  products  are  required.  Furthermore,  expressions  for  the  internal 
energy  are  also  needed  to  test  the  conservation  of  energy  relation  (see  equation  2). 
These  expressions  are  dependent  upon  the  equation  of  state  chosen  to  represent  the 
detonation  products.  Mader  has  derived  the  appropriate  expressions  for  the  BKW 
equation  of  state. 

The  BKW  EOS  has  the  following  form  [1]; 

PV 

— i-  =  1  +  xepx  =  F(x)  (10) 

RT 

Kk 

where  x  = -  (11) 

Vg(T  +  6)a 

N  f 

and  k  =  Xx-kj  (12) 

i=l 


f  # 

In  these  expressions  Vg  is  the  molar  volume  of  the  gaseous  products,  Xj  is  the  mole 

fraction  of  the  ith  component  of  the  equilibrium  composition,  and  kj  is  the  co-volume 
of  the  ith  gaseous  component  (which  is  effectively  an  estimate  of  the  molecular 
volume).  a,p,K  and  0  are  equation  of  state  parameters  whose  values  depend  on  which 
particular  model  is  adopted.  Note  that  the  BKW  equation  of  state  only  applies  to  the 
gaseous  detonation  products.  For  solid  products,  such  as  carbon  or  AI2O3,  we  use 
the  Cowan-Fickett  equation  of  state,  which  has  the  form  [16]; 

P  -  PiCn) +  a(T|)T  +  b(r[)T2  (13) 

where  T)  =  p  /  p0  is  the  compression  of  the  solid  material  relative  to  its  normal  crystal 
density,  and  p:  a  and  b  are  polynomial  expressions  in  r\. 

* 

For  gaseous  detonation  products  described  by  the  BKW  EOS  Fj  is  given  by 


Fi  = 


(  F°-H ^ 


RT 


Ji 


W)i 

RT 


+  Ini 


x  P 


pfr-l  -  k-  - 

+^-ir^-\nF(X)+^-(F(x)-l) 
p  k 


(14) 


In  equation  (14),  (F°  -  H®);  is  the  free  energy  of  species  i  if  i  is  regarded  as  a  perfect 
gas  and  N  is  the  total  number  of  gaseous  products.  The  free  energy  term  is  calculated 
using  the  following  expression: 
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F°  -  H°n 


b.T  cgT 2  dj3  egT 
va*+2+3  +  4  5y 


ric„ 


+  - 


(15) 


where  ag/  bg/  cg/  dg/  eg  and  ricg  are  constants. 

The  expression  for  the  free  energy  of  the  solid  components  is 


F*- 


^F°  -H°A 


RT 


A 


I  (Ho)s  ,  Fs 
RT  RT 


(16) 


where 


Fs=  Mr 


PVS  -  P0  v0 


asV  +  bs/nV 


's  ds 


V  2V2  3V3 


+  (AjV  +  A2  In  V)Ty  + 


QV  + 


c3v3> 

™  2 

V0 

+  — — 

Ty 

3 

7 

- 

Vo, 

(17) 


where  Tv  is  the  temperature  in  volts,  Mr  is  the  molecular  mass  in  g/mol,  V0  is  the 
S.T.P.  volume  of  the  solid  in  cm3/g  and  Vs  is  the  solid  volume  in  cm3/ g  at  the  relevant 
P  and  T.  Ah  A2,  C2,  C2,  C3,  as,  bs,  cs,  ds  and  es  are  the  Cowan  equation  of  state 

constants  specific  for  each  solid  product.  ((F°  -H0  )/T)s  is  computed  as  per  equation 
15. 


The  internal  energy  is  given  by  [1]: 


-  RT  + 


ocRT2 
T  +  0 


(F(x)  -  1) 


Esi  =  (H°  -  Hg)i  +  (AHf  )°+REg 


H°  -  H°  = 


bgT^ 


+ 


2cgT- 


+ 


3d„T4  4egT5 

— ± H - 5 —  +  nC 


(18) 


(19) 


(20) 


where  H°  -  H°  is  the  enthalpy  of  the  product  in  the  standard  state  at  T  Kelvin  minus 
the  enthalpy  of  the  product  in  the  standard  state  at  0  Kelvin. 


Es  =  Mr 


QV  + 


C2v2  t  C3V 


-  + 


Ty- 


( 


ae V  +  b  AnV  — -  - 


V 


V  2V2 


(21) 


and  Ftotal -  XxjEg(i)  (22) 

solids 
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2.4  The  AMRL  Code 

Given  the  ability  to: 

(i)  compute  the  equilibrium  composition  at  a  given  P  and  T  by  minimizing  the 
free  energy. 

(ii)  compute  the  internal  energy  of  this  equilibrium  mixture  at  P  and  T 

then  the  Hugoniot  equation  and  the  CJ  condition  can  be  solved  to  yield  the  CJ  point. 

Use  of  the  CJ  condition  in  the  form  given  by  equation  (1)  is  inconvenient  when 
coupled  with  the  form  of  equations  (18)  through  (22)  giving  internal  energy.  By 
recalling  that 


fsp)  . 

(dp] 

'oT 

II 

Is  V) 

at  the  CJ  point,  we  can  use  the  CJ  condition  in  the  form 


(  p  -p\ 
0 

U  v)r 

V  -V 
\vo  v  y 

(23) 


(24) 


The  CJ  point  is  then  determined  via  the  following  sequence: 

(i)  make  an  initial  estimate  of  P,  T,  and  the  equilibrium  composition. 

(ii)  minimize  the  free  energy  to  obtain  the  correct  composition  for  the  given  P 
and  T. 

(iii)  at  this  stage  the  variables  P,  T  do  not  necessarily  lie  on  the  Hugoniot,  so 
equation  (2)  is  solved  to  obtain  a  new  value  of  T.  The  composition  at  this 
new  value  of  T  is  incorrect  now,  and  so  the  free  energy  must  be  minimized 
again  to  obtain  the  correct  composition  at  the  new  value  of  T. 

(iv)  after  the  above  iteration  we  have  a  value  of  P  and  T  which  lie  on  the 
Hugoniot  and  for  which  the  composition  is  valid,  but  the  variables  do  not 
necessarily  satisfy  the  CJ  condition.  Hence  equation  (24)  is  now  solved  to 
determine  a  new  value  of  P  which  satisfies  the  CJ  condition.  The 
equilibrium  composition  is  now  no  longer  valid  at  the  new  value  of  P,  and 
so  the  sequence,  (ii),  (iii)  and  (iv)  is  repeated. 


As  can  be  seen  from  the  above  procedure,  the  calculation  of  the  CJ  state  is  an  iterative 
process,  and  the  iteration  is  continued  until  the  difference  between  the  current  value  of 
the  variables  and  the  previous  values  falls  below  some  predetermined  value. 

Once  the  CJ  state  has  been  determined,  the  detonation  velocity  at  the  CJ  point  can  be 
calculated  using: 
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Da  =  CfV0 


Pcj  -  Pq 
Vo  -  VCJ 


(25) 


where  Pq  and  Vq  are  the  pressure  in  Mbar  and  specific  volume  in  cm3/g 
respectively.  Cf  is  a  conversion  factor  (Cpl.OxlO4)  to  give  the  detonation  velocity  in 
m/s.  A  computer  program  was  written  (called  SDA»FOR)  to  implement  the  above 
equations  and  was  used  to  calculate  the  detonation  properties  of  a  wide  range  of 
explosives  and  explosive  compositions. 


2.5  Code  Validation 

Table  1  contains  a  listing  of  the  explosives  and  the  calculated  CJ  states,  as  well  as  a 
comparison  with  the  results  obtained  from  Mader's  BKW  code.  For  these  calculations 
a  single  set  of  BKW  equation  of  state  parameters  was  used.  These  values  are  a=0.5, 
(3=0.16,  k=10.90978  and  6=400;  this  set  is  referred  to  as  the  RDX  parameter  set  by 
Mader  [1].  The  agreement  between  SDA»FOR  and  BKW  is  excellent  for  all  explosives 
and  explosive  compositions  shown  in  Table  1.  Typically  the  percentage  difference 
between  the  BKW  and  SDA  calculated  detonation  velocities  is  less  than  0.5%.  For  the 
detonation  pressure  a  percentage  difference  of  5%  is  obtained  for  PETN  at  a  loading 
density  of  0.5  g/cm3  and  at  a  loading  density  of  1.0  g/cm3  the  same  explosive  has  a 
percentage  difference  of  3%.  Apart  from  these  two  cases  the  percentage  difference  in 
the  detonation  pressure  is  less  than  1%.  Suceska  [21]  has  written  a  code  (EXPL05)  to 
perform  calculations  of  detonation  parameters  using  a  BKW  equation  of  state  and  has 
compared  results  of  this  code  with  results  obtained  from  the  Mader  BKW  program.  In 
this  comparison  Suceska  also  observes  small  differences  between  the  two  programs 
and  suggests  that  the  differences  could  be  due  to  the  use  of  different  values  for  the 
standard  thermodynamic  functions  or  due  to  a  difference  in  the  way  EXPL05 
computes  solid  product  thermodynamic  functions.  The  minor  variations  observed 
between  the  SDA  code  and  the  Mader  BKW  code  can  be  attributed  to  coding 
differences  between  the  two  programs.  For  example,  in  the  SDA  code  an  iterative 
technique  is  used  to  determine  the  CJ  point  whereas  in  the  Mader  BKW  code  the  CJ 
point  is  determined  by  approximating  the  detonation  velocity  as  a  parabolic  function 
in  P  and  then  finding  the  minimum  of  this  parabola  to  give  the  minimum  detonation 
velocity  and  hence  the  CJ  state.  This  and  other  differences  in  numerical  techniques  can 
explain  the  minor  variations  noted  above. 


8 


DSTO-TR-0226 


Table  1:  Comparison  of  computed  C]  parameters  for  selected  explosives  and 
explosives  compositions.  The  values  in  the  BKW  columns  were  computed  by  the 
Mader  BKW  code  (see  Mader  [1])  and  the  values  in  the  columns  headed  by  SDA  were 
computed  using  the  SDA.FOR  code  written  for  this  work. 


EXPLOSIVE 

Pcj 

(Mbar) 

Tcj 

(K) 

Dcj 

(m/s) 

BKW 

SDA 

BKW 

SDA 

BKW 

SDA 

RDX  p  =  1.8 

0.347 

0.345 

2587 

2588 

8754 

8711 

TNT  p  =  1.64 

0.213 

0.213 

2829 

2825 

7197 

7168 

HMX  p  =  1.9 

0.395 

0.394 

2364 

2364 

9159 

9121 

PETN  p  =  1.67 

0.280 

0.279 

3018 

3014 

8056 

8024 

p  =  1.0 

0.101 

0.0978 

3970 

3958 

5947 

5929 

p  =  0.5 

0.0303 

0.0286 

4493 

4599 

4313 

4308 

TATB  p=  1.895 

0.326 

0.325 

1887 

1890 

8411 

8365 

PADP  p  =  1.86 

0.300 

0.298 

3112 

3112 

7971 

7931 

HNS  p  =1.74 

0.241 

0.241 

3059 

3057 

7410 

7377 

RDX/TNT  64/36  (Comp  B) 

p  =  1.713 

0.284 

0.282 

2763 

2764 

8084 

8037 

RDX/TNT/Wax  48.9/46.1/5 

p  =  1.62 

0.237 

0.236 

2741 

2738 

7609 

7576 

TNT/FETN  50/50 

P*l-65 _ 

0.257 

0.256 

3239 

3235 

7740 

7707 

RDX=cyclotrimethylene  trinitramine,  TNT=2, 4/6-trinitrotoluene,  HMX=cyclotetramethylene 
tetranitramine,  PETN=pentaerythritol  tetranitrate,  TATB=l,3,5-triamino-2, 4,6- trinitrobenzene, 
PADP=2,6-bis(picrylazo)-3,5-dinitropyridine. 
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3.  The  Probabilistic  Algorithm  for  Energy 
Minimization 


One  of  the  disadvantages  of  the  method  of  steepest  descent  described  in  the  previous 
section  is  the  need  to  provide  an  initial  estimate  of  the  equilibrium  composition  which, 
in  some  cases,  needs  to  be  reasonably  close  to  the  "correct"  equilibrium  composition. 
A  further  problem  is  the  need  to  ensure  that  all  the  Xj  remain  positive  for  each 
iteration.  With  the  method  of  steepest  descent  it  was  found  that  if  any  of  the  assumed 
products  was  particularly  unfavourable  (ie.  had  a  very  large  free  energy)  then  the 
method  would  try  to  make  the  Xj  of  that  product  negative.  This  can  be  avoided  by 
employing  a  scaling  factor,  but  in  some  cases  the  scaling  factor  is  so  small  that  the 
composition  barely  changes  from  one  cycle  to  the  next  and  the  true  equilibrium 
composition  cannot  be  obtained. 

White  et  al.  [17]  described  an  alternative  method  for  minimizing  the  free  energy 
which  is  based  on  a  linear  programming  technique  and  automatically  ensures  that  all 
xj  remain  positive  during  the  course  of  the  calculation.  This  method  was  considered 
for  use  in  the  SDA«FOR  code  but  was  found  to  be  unsuitable  because  in  the 
formulation  of  the  problem  it  is  necessary  for  the  total  free  energy  function  to  be 
expressed  as  a  linear  function  of  the  product  composition  vector,  and  with  BKW  (and 
other  complex  EOS)  this  is  not  the  case. 


3.1  Description  of  Algorithm 

An  alternative  approach  to  minimizing  functions  is  to  use  a  probabilistic  or  stochastic 
methodology.  In  this  general  approach,  a  composition  with  Xj  >  0  is  guessed  and  the 
free  energy  is  computed.  This  process  is  continued  until  a  pre-determined  number  of 
guesses  has  been  achieved.  At  this  stage  the  equilibrium  composition  is  that 
composition  amongst  all  the  guesses  that  yielded  the  lowest  free  energy.  More 
sophisticated  versions  employ  a  weighting  function  approach  to  improve 
convergence.  The  advantages  of  a  stochastic  method  include  the  guarantee  that  xj  >  0, 
and  that  a  global  maximum/ minimum  will  be  found  (deterministic  algorithms  such  as 
the  method  of  steepest  descent  find  local  maxima /minima  and  do  not  guarantee 
finding  a  global  one).  A  potential  disadvantage  of  a  probabilistic  approach  is  that  the 
total  computing  time  will  be  longer. 

A  stochastic  method  presented  by  Benke  &  Skinner  [18]  uses  an  'adaptive  probabilistic 
algorithm'  for  locating  global  optima  of  multivariate  functions.  Let  f  be  a  function  of 
some  vector  X  [in  the  application  to  an  equilibrium  code  f  =  free  energy  and  X  = 
product  composition  vector }.  The  algorithm  can  then  be  represented  by  the  following 
pseudo  code  {virtually  identical  to  that  shown  by  Benke  &  Skinner). 
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1.  Generate  Xj  randomly  within  constraints 

2.  f]  =  function  (X|) 

3.  FOR  number  of  guesses  DO 

generate  X2  randomly  within  constraints 
f2  =  function  (X2) 

X3  =  (f|  X|  +  f2  X2)  /  (fj  +  f2) 
f3  =  function  (X3) 

IF  fj  =  min  (fj,  f2,  f3)  THEN 
do  nothing 

ELSE 

IF  f2  =  min  (fj,  f2,  f3)  THEN 
X1  =  X2 

h  =  f2 

ELSE 

Xi  =  X3 
h  =  f3 

ENDIF 

ENDIF 

ENDFOR 

4.  X|  is  optimum  parameter  vector 
fj  is  optimum  value 

5.  END  algorithm 

The  weighted  mean  X3  is  the  adaptive  part  of  the  algorithm.  Xx  is  the  current  best 
guess  at  any  stage  of  the  algorithm. 


3.2  Implementation  in  AMRL  Code 

In  order  to  apply  this  algorithm  to  the  problem  of  determining  equilibrium 
compositions,  we  need  a  method  of  generating  random  composition  vectors  X  that 
sample  the  allowed  solution  space  of  X  evenly.  The  first  attempt  was  based  on  the 
premise  that  allowed  composition  vectors  can  be  obtained  by  taking  linear 
combinations  of  basis  vectors  that  span  the  allowed  solution  space.  Unfortunately, 
although  such  basis  vectors  could  be  determined,  it  was  found  that  linear  combination 
of  these  vectors  gave  rise  to  x,  <  0. 

Failing  this,  a  new  approach  was  implemented.  Consider  a  list  of  n  products  X  =  (xj, 
x2,....xn)  where  each  product  contains  one  or  more  of  m  elements.  Let  A  be  the 
product  composition  matrix  where  ajj  specifies  the  number  of  atoms  of  type  j  which 
are  in  product  i.  Let  B  be  the  initial  or  reactant  vector  where  bj  specifies  the  initial 
amount  of  atom  type  j.  Choose  a  product,  number  P,  at  random  from  the  list  and 
assign  it  a  random  composition,  x~,  between  0  and  an  upper  limit  given  by  min 
(bj/apj)  thus  ensuring  that  none  of  the  bj  will  go  negative.  Update  the  bj  values. 
Now  select  another  product  at  random  until  all  products  have  been  processed.  It  is 
possible  that  at  the  end  of  the  process  not  all  of  the  bj  =  0;  in  this  case  the  composition 
determined  is  invalid  and  is  thus  rejected.  Another  feature  is  added  to  the  algorithm 
to  reduce  the  number  of  rejected  compositions.  When  a  product  is  selected  a  check  is 
made  to  see  if  it  is  the  last  product  left  that  contains  a  particular  element.  If  so  then 
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the  product  is  assigned  a  composition  such  that  all  of  the  element  in  question  is  used 
up.  The  algorithm  can  be  expressed  as: 

1.  Set  up  list  of  products 

2.  Set  up  counters  indicating  number  of  products  for  each  element 
ie  count  (j)  =  number  of  products  containing  j'th  element. 

3.  WHILE  (more  products  in  list)  DO 

Select  product  at  random. 

Determine  maximum  x  for  this  product: 
xmax  =  min  (tot(j)/prod(j);  prod  (j)  *  0) 

IF  (not  the  last  product  containing  any  element)  then 
x  =  random  number  between  0.0  and  xmax 

ELSE 

x  =  xmax 

ENDIF 

IF  (x  <  0)  THEN  start  again  {go  to  1} 

FOR  (each  element)  DO 

IF  (product  contains  element;  ie  prod(j)  *  0)  THEN 
tot(j)  =  tot(j)  -  x  *  prod  (j) 
count  (j)  =  count  (j)  - 1 

ENDIF 

ENDFOR 

Remove  product  from  the  list 
ENDWHILE 

4.  IF  (any  tot(j)  *0)  THEN  start  again  {go  to  1} 

5.  Finished. 

Note:  prod(j)  =  number  of  atoms  of  type  j  in  the  chosen  product 

tot(j)  s  number  of  atoms  of  type  j  left;  initially  the  tot(j)  value  reflect  the 
reactant  composition. 

A  new  program  PEA«FOR  (for  Probabilistic  Equilibrium  Algorithm)  was  written  to 
incorporate  the  minimization  method  of  Benke  &  Skinner.  The  program  uses  the  BKW 
EOS,  a  bisection  method  to  solve  the  Hugoniot  equation,  and  a  VMS  random  number 
generator. 


3.3  Code  Validation 

Table  2  shows  a  comparison  of  the  results  obtained  with  PEA*FOR  and  with  the  BKW 
code.  Once  again  a  single  set  of  BKW  parameters,  the  RDX  parameters  cited  earlier,  is 
used  for  the  values  quoted  in  Table  2.  These  results  indicate  that  the  program  can 
determine  CJ  values  for  explosives  that  are  in  close  accord  with  the  results  of  Mader's 
BKW  code.  It  is  also  evident  that  the  code  is  much  more  time  intensive  (all  runs  were 
performed  on  a  VAX  8700),  but  the  location  of  a  global  minimum  of  free  energy 
without  a  good  initial  guess  is  virtually  guaranteed.  In  addition,  no  problems  with  Xj 
<  0  arise  with  this  method. 
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Table  2:  Comparison  of  computed  CJ  parameters  calculated  using  the  conventional 
method  of  steepest  descent  (BKIN  column)  and  a  probabilistic  algorithm  for 
determining  equilibrium  compositions. 


Probabilistic  Code 

BKW 

Explosive 

Pq 

(Mbar) 

Dq 

(m/s) 

Tq 

(K) 

CPU 

time 

#  of 
guesses 

Pq 

(Mbar) 

Dq 

(m/s) 

Tq 

(K) 

Comp  B 
p  =  1.713 

0.289 

8081 

2779 

58  hr 

29  min 

100,000 

0.284 

8084 

2763 

PETN 
p  =  1.67 

0.285 

8076 

2898 

2  hr 

31  min 

5,000 

0.280 

8056 

3018 

0.285 

8078 

2957 

5  hr 

2  min 

10,000 

0.286 

8081 

3015 

49  hr 

55  min 

100,000 

HMX 
p  =  1.90 

0.402 

9197 

2274 

3  hr 

11  min 

5,000 

0.395 

9195 

2364 

0.404 

9207 

2330 

6  hr 

5  min 

10,000 

0.402 

9183 

2366 

61  hr 

21  min 

100,000 

TNT 
p  =  1.64 

0.246 

7317 

2949 

2  hr 

44  min 

5,000 

0.213 

7197 

2829 

0.216 

7191 

2866 

4  hr 

42  min 

10,000 

0.218 

7222 

2831 

47  hr 

13  min 

100,000 

RDX 

p  =  1.8 

0.353 

8790 

2510 

2  hr 

23  min 

5,000 

0.347 

8754 

2587 

0.355 

8853 

2573 

4  hr 

42  min 

10,000 

0.353 

8836 

2588 

51  hr 

17  min 

100,000 
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4.  Application  to  PBXN-111 


4.1  The  Method  of  Steepest  Descent 

Prior  to  computing  detonation  parameters  of  the  non-ideal  explosive  formulation 
PBXN-111,  a  description  of  the  formulation  is  required.  In  particular  the  equilibrium 
code  SDA.FOR  requires  the  elemental  breakdown,  heat  of  formation  and  density  of 
the  formulation.  PBXN-111  comprises  20%  RDX  (cyclotrimethylene  trinitramine),  43% 
AP  (ammonium  perchlorate),  25%  A1  (aluminium)  and  12%  of  a  HTPB  based 
polyurethane  binder.  The  binder  can  be  broken  down  into  its  components;  47.7% 
HTPB  (hydroxy  terminated  polybutadiene),  47.7%  IDP  (isodecyl  pelargonate)  and 
4.6%  IPDI  (isophorone  diisocyanate).  The  heat  of  formation  of  the  formulation  is 
calculated  by  summing  the  contribution  of  each  of  the  components  as  shown  in  Table 
3.  The  elemental  composition  is  obtained  in  a  similar  fashion  by  summing  the 
contribution  of  each  component  to  the  amount  of  each  element  as  shown  in  Table  4. 
Given  an  explosive  loading  density  of  p=1.80  g/cm3  for  PBXN-111,  the  CJ  parameters 
can  be  determined  using  the  SDA  code  and  an  appropriate  set  of  products. 


Table  3:  Nominal  formulation  of  PBXN-111  and  relevant  heats  of  formation. 


Component 

AHf(0  K) 
kcal/mol 

Molecular 

Formula 

Mr 

g/mol 

Weight 

Fraction 

Contribution 
to  E0 
cal/mol 

RDX 

33.97 

CaHfiNfiOfi 

222.11 

0.20 

3058.8 

AP 

-70.73 

nh4cio4 

117.489 

0.43 

-25886.6 

Al 

0 

Al 

26.98 

0.25 

0 

HTPB 

-0.038* 

c7,33h11°0.083 

100.0 

0.0572 

-2.2 

IDP 

212.8* 

c19h38°2 

298.51 

0.0572 

4077.6 

IPDI 

-88.8* 

C12H18N2°2 

222.29 

0.0056 

-223.7 

Total 

-18976.1 

*  heats  of  formation  for  these  components  are  not  available  at  0  K  soAHf  at  298  K  is  used  instead.  This  is 
not  expected  to  significantly  alter  the  results. 


The  parameter  set  used  in  the  BKW  equation  of  state  was  obtained  by  adjusting  the 
parameters  to  fit  observed  CJ  values  for  a  number  of  explosives.  This  procedure  has 
been  performed  a  number  of  times  by  different  workers  and  thus  different  parameter 
sets  are  available.  In  Table  5  the  calculated  detonation  properties  of  PBXN-111  are 
shown  using  SDA.FOR  and  different  equation  of  state  parameter  sets.  Although  there 
are  differences  between  the  calculated  detonation  velocity  and  detonation  pressure  for 
the  four  parameter  sets  shown  in  Table  5,  the  calculated  values  are  all  significantly 
higher  than  those  obtained  from  experimental  observations.  Forbes  [5]  has  reported  a 
CJ  pressure  for  PBXN-111  of  12.2  GPa  (0.120  Mbar),  which  is  less  than  half  the 
pressure  calculated  using  SDA.FOR  and  any  of  the  four  parameter  sets.  The  infinite 
diameter  detonation  velocity  obtained  by  Forbes  using  a  linear  extrapolation  technique 
is  6195  m/s,  while  Bocksteiner  et  al.  [8],  using  the  same  linear  extrapolation,  report  a 
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lower  value  of  5913  m/s.  Both  experimental  estimates  are  very  much  lower  than  the 
values  shown  in  Table  5.  Bocksteiner  et.  al.  have  shown  that  an  elliptical  fit  of  the 
detonation  velocity  data  is  superior  to  a  linear  fit.  Using  an  elliptical  extrapolation 
they  obtain  an  infinite  diameter  detonation  velocity  of  5641  m/ s  for  Australian  PBXN- 
111  and  a  value  of  5760  m/s  for  US  PBXN-111. 


Table  4:  Elemental  composition  of  PBXN-111. 


Component 

Atom 

C 

H 

N 

O 

Cl 

Al 

RDX 

0.270 

0.540 

0.540 

0.540 

- 

- 

AP 

- 

1.464 

0.366 

1.464 

0.366 

- 

Al 

- 

- 

- 

- 

- 

0.927 

HTPB 

0.419 

0.629 

- 

0.005 

- 

- 

IDP 

0.364 

0.728 

- 

0.038 

- 

- 

IPDI 

0.030 

0.045 

0.005 

0.005 

- 

- 

Total 

1.083 

3.406 

0.911 

2.052 

0.366 

0.927 

Two  possibilities  for  the  large  discrepancy  between  calculated  and  experimental 
estimates  of  the  CJ  pressure  and  detonation  velocity  will  be  considered  here.  Firstly,  it 
is  well  known  that  certain  explosive  formulations  behave  non-ideally  [23],  and  thus 
calculations  based  upon  the  steady  state  theory  developed  by  Chapmen  and  Jouguet 
cannot  reproduce  experimental  values.  Secondly,  the  parameters  in  the  BKW  equation 
of  state  are  determined  by  fitting  calculations  to  experimental  results.  Consequently 
the  equation  of  state  is  "customized"  for  the  particular  explosives  used  in  the  fitting 
process.  This  being  the  case,  one  would  expect  to  obtain  reasonable  results  for 
explosives  similar  to  the  explosives  used  in  the  fit.  However,  if  the  explosive  under 
consideration  is  very  different  then  the  reliability  of  the  calculated  results  is  reduced. 
A  more  general,  non-parameterized  equation  of  state  would  allow  greater  flexibility  in 
predicting  detonation  parameters  of  new  explosives/ formulations. 


Table  5:  Calculated  CJ  state  of  PBXN-111. 


RDX  parameters 

TNT  parameters 

BKWR  set  [23] 

Hobbs  and  Baer  [22] 

T  (K) 

5229 

5291 

5494 

5823 

P  (Mbar) 

0.294 

0.299 

0.305 

0.289 

D  (m/s) 

8186 

8090 

8337 

7909 

a 

0.5 

0.5 

0.5 

0.5 

p 

0.16 

0.09585 

0.176 

0.174 

e 

400 

400 

1850 

5160 

K 

10.90978 

12.685 

11.8 

11.85 
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Mader  [1]  defines  a  non-ideal  explosive  as  having  a  CJ  pressure,  velocity  or  expansion 
isentrope  significantly  different  from  those  expected  from  equilibrium  steady  state 
calculations.  PBXN-111  certainly  satisfies  this  definition.  There  are  a  number  of 
reasons  why  an  explosive  may  behave  non-ideally.  For  example  the  reaction  of  one  or 
more  components  of  the  explosive  may  occur  on  a  timescale  too  long  to  support  the 
propagation  of  the  shock  front,  the  explosive  geometry  may  lead  to  rarefactions  that 
reduce  the  energy  transferred  to  the  shock  front  or  the  magnitude  of  the  initiating 
pulse  can  affect  the  detonation  state  achieved  by  a  non-ideal  explosive.  Using  the 
chemical  equilibrium  code  BKW,  Johnson,  Mader  and  Goldstein  [10]  have  examined 
the  influence  of  partial  reaction  of  one  of  the  reactants  on  the  CJ  parameters  of  Amatex 
40.  Their  results  showed  that  by  assuming  50%  of  the  AN  (ammonium  nitrate) 
remains  inert  then  the  BKW  calculation  closely  matches  the  experimental  value. 

Using  this  idea,  a  series  of  calculations  where  the  extent  of  reaction  of  one  or  more  of 
the  reactants  is  limited  was  undertaken  with  the  SDA  code.  The  aim  of  these 
calculations  was  to  ascertain  what  degree  of  reaction  of  the  A1  and  AP  would 
reproduce  experimental  detonation  velocity  and  pressure  for  PBXN-111.  The  first  set 
of  calculations  involved  taking  RDX  and  successively  adding  an  inert  diluent  to  gauge 
the  effect  of  this  inert  diluent  on  the  RDX  explosive  performance.  In  these  calculations 
A1  was  added  as  the  inert  diluent.  Although  A1  is  reactive,  for  the  purpose  of  this 
exercise  it  was  assumed  to  be  inert.  Figure  I  shows  that  both  the  detonation  velocity 
and  pressure  decrease  with  increasing  percentage  of  the  A1  diluent.  Thus  the  addition 
of  a  component  that  is  either  unreactive  or  does  not  react  on  a  timescale  appropriate  for 
sustaining  a  detonation  will  reduce  the  explosive  performance  as  calculated  by  the 
equilibrium  code. 

Similar  calculations  were  performed  for  a  PBXN-111  like  formulation.  For  these 
calculations  the  formulation  was  simplified  by  removing  the  binder,  giving  a  mixture 
containing  20%  RDX,  49%  AP  and  31%  Al.  In  these  calculations,  the  A1  was  assumed 
to  be  inert  and  the  amount  of  AP  that  participated  in  the  reaction  was  varied.  The  inert 
AP  was  assumed  to  behave  similarly  to  inert  AN  (ammonium  nitrate).  This 
assumption  was  necessary  since  appropriate  thermodynamic  and  solid  equation  of 
state  parameters  were  not  available  for  AP  but  were  readily  available  for  AN.  The 
results  of  these  calculations  are  shown  in  Figure  2.  Once  again  the  calculations  reveal 
that  both  the  detonation  velocity  and  detonation  pressure  decrease  as  the  amount  of 
inert  material  is  increased.  The  result  obtained  when  all  of  the  AP  is  assumed  to  be 
inert  gives  DCj=6475  m/s  and  PCj=0.143  Mbar.  This  result  is  still  somewhat  higher 
than  the  values  obtained  by  Forbes  [5];  Dq=6195  m/s  and  Pq=0.120  Mbar.  Based  on 
the  comparison  of  the  experimental  results  and  the  results  of  the  calculations  presented 
in  Figure  2,  it  appears  that  RDX  is  the  major  driving  force  in  the  detonation  of  PBXN- 
111  and  that  AP  and  Al  provide  very  little,  if  any,  contribution  to  the  propagation  of 
the  shock  wave. 

4.2  The  Probabilistic  Algorithm 

The  previous  section  has  shown  that  attempts  to  predict  the  detonation  velocity  and 
pressure  of  PBXN-111  using  the  BKW  EOS  and  the  method  of  steepest  descent 
minimization  algorithm  lead  to  values  which  are  significantly  higher  than  the 
experimental  values.  One  possible  cause  of  this  discrepancy,  which  was  discussed  in 
the  previous  section,  is  that  the  BKW  EOS  has  been  parametrised  on  a  data  base 
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containing  primarily  ideal  explosives,  and  hence  is  unsuited  for  use  with  highly  non¬ 
ideal  explosives.  The  previous  section  showed  that  restricting  the  extent  of  reaction  of 
both  the  AP  and  A1  considerably  reduced  both  the  detonation  velocity  and  pressure, 
but  the  best  estimates  obtained  were  still  significantly  higher  than  the  experimental 
results.  If  the  minimization  algorithm  is  finding  the  global  minimum  in  the  free 
energy,  then  the  above  results  suggest  that  use  of  a  more  appropriate  equation  of  state 
is  required.  We  consider  this  possibility  in  the  next  section. 


Another  possible  cause  of  the  discrepancy,  which  was  mentioned  in  the  Introduction, 
is  that  the  algorithm  which  finds  the  minimum  in  the  free  energy  is  failing  to  find  the 
true  free  energy  minimum.  To  investigate  this  possibility  we  now  use  the  probabalistic 
algorithm  in  the  AMRL  code  to  calculate  the  CJ  state  of  PBXN-111.  Unfortunately  the 
amount  of  computer  time  required  to  perform  this  calculation  is  considerably  higher 
than  the  time  needed  for  the  simpler  CHNO  explosives  .  PBXN-111  contains  two  extra 
elements,  namely  Cl  and  Al,  in  comparison  to  standard  CHNO  explosives. 
Consequently  there  are  more  products  in  the  composition  and  this  considerably 
increases  the  computer  time  required  to  find  the  equilibrium  composition.  Table  6 
summarise  the  results  obtained  when  the  PEA-FOR  program  was  used  to  predict  the 
CJ  state  of  PBXN-111  using  the  RDX  set  of  parameters.  The  composition  and  CJ 
parameters  are  shown  for  a  number  of  initial  guesses,  varying  from  1,000  to  100,000. 
Comparing  the  results  for  100,000  guesses  with  the  CJ  parameters  calculated  using  the 
method  of  steepest  descent  and  the  RDX  parameters  shows  an  approximate  2% 
difference  in  CJ  pressure,  4%  difference  in  CJ  velocity,  and  1.5%  difference  in  CJ 
temperature.  This  indicates  that  the  probabalistic  algorithm  is  finding  a  different 
minimum  to  the  one  found  by  the  method  of  steepest  descent.  However,  the 
differences  between  the  results  predicted  by  the  two  algorithms  are  small  compared  to 
the  differences  between  the  experimental  values  and  the  results  calculated  by  either  of 
the  minimization  algorithms,  and  so  we  conclude  the  BKW  EOS  is  inappropriate  for 
the  description  of  the  CJ  state  of  explosives  such  as  PBXN-111. 

Table  6:  Comparison  of  computed  C]  parameters  and  compositions  calculated  for 
PBXN-111  with  different  number  of  guesses  using  the  probabilistic  algorithm  for 
determining  equilibrium  compositions. 


Number  of  guesses 

Composition 

1000 

2000 

3000 

10,000 

100,000 

HzO 

0.6435 

0.6435 

0.6480 

0.6432 

0.6481 

CCh 

0.0023 

0.0023 

0.0011 

0.0013 

0.0004 

Nz 

0.1606 

0.1606 

0.1603 

0.1600 

0.1595 

NHa 

0.5898 

0.5898 

0.5904 

0.5909 

0.5920 

CO 

0.0135 

0.0135 

0.0110 

0.0157 

0.0126 

HC1 

0.3497 

0.3497 

0.3622 

0.3467 

0.3337 

Ch 

0.0081 

0.0081 

0.0139 

0.0097 

0.0161 

Cm 

1.0673 

1.0673 

1.0709 

1.0660 

1.0700 

AI2O3M 

0.4635 

0.4635 

0.4635 

0.4635 

0.4635 

Cl  parameters 

Pci  (Mbar) 

0.257 

0.257 

0.290 

0.258 

0.300 

Tci  (K) 

5202 

5202 

5286 

5203 

5309 

Dcj  (m/s) 

7822 

7822 

7861 

7821 

7884 

CPU  time 

2hrs  44  min1 

8hrs  3  min2 

9hrs  37  min2 

39hrs 

llmin2 

322hrs 

48min2 

1:  Run  on  a  HP  715  /75  2:  Run  on  a  HP  715/50 
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5.  The  Effect  of  Intermolecular  Equations  of  State 


5.1  The  JCZ3  Equation  of  State 

As  discussed  in  the  Introduction,  previous  work  has  shown  that  highly  non-ideal 
composite  explosives  are  best  described  using  an  EOS  based  on  an  intermolecular 
potential.  The  JCZ3  equation  of  state,  as  described  by  Cowperthwaite  and  Zwisler  [19], 
satisfies  this  criterion,  and  has  been  shown  to  be  applicable  to  several  non-ideal 
explosives  [13,14].  The  advantage  of  JCZ3  over  BKW,  for  example,  is  that  it  does  not 
contain  any  parameters  which  must  be  determined  from  an  empirical  fit.  JCZ3  would 
therefore  be  expected  to  be  equally  applicable  to  a  wide  range  of  explosives,  whereas 
BKW  is  only  useful  for  those  classes  of  explosives  used  in  the  fitting  process.  Hence,  in 
this  section,  we  describe  the  implementation  of  JCZ3  in  the  SDA*FOR  code 


5.2  Implementation  of  JCZ3  EOS  in  the  AMRL  Code 

Incorporating  the  JCZ3  equation  of  state  into  the  SDA  code  involves  deriving 
expressions  for  the  free  energy  and  internal  energy  of  a  mixture  of  gaseous  products. 
For  a  mixture  of  n  moles  of  s  fluid  species: 

nRT 

P  =  p0(V)  +  G(V,T)—  (26) 

where  p0(V)  =  lattice  pressure  along  the  zero  degree  isotherm,  G(V,T)  is  the  factor  that 
accounts  for  the  thermal  contribution  to  the  pressure  from  intermolecular  forces,  and  R 
is  the  universal  gas  constant. 

po(V)  has  the  form  [19] 


p  (V)  =  -^o 

dV 

where  Eo(V)  is  the  volume  potential  of  a  face-centered  cubic  lattice,  given  by: 


E0<v)  =  ^q 


B, 


exp 


/-/ 


vvy 


B 


m 


m 

r„*xj 


m 


where 


(27) 


e0  =  -  E  I  nin  jejj 


n  i  j 


(28) 
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V*  =-EIniIljV*j 
n  i  j 


* 


rU  = 


fj  +rj 
2 


(29) 

(30) 

(31) 

(32) 


and  B /  =  13.99166,  Bm  =  14.45392,  /  =  repulsive  exponent  =  13.5,  m  =  attractive 
exponent  =  6,  &  =  potential  well  depth,  n  =  equilibrium  distance,  and  N  =  Avogadro's 
number.  The  G  factor  is  given  by: 


where 


G(V,T)  =  1-y 


f  =  fg+fs 


(33) 


(34) 


fs  ^  nRT  (/  -  mj 


m  }l_ 

71 


if 


V' 


-  2  exp 


l- l 


/V\  3 


VV2 


fg  =  l  +  a!y  +  a2y2  +a3y3 


(35) 

(36) 


y 


(37) 


F  =  C! 


In 


T(/-m) 
m(e0  /  nR) 


(38) 


Cj  =  c  +  / 


(39) 


and  c  s  Euler's  constant  =  0.57722.  The  chemical  potential  of  the  i'th  species  can  be 
obtained  using: 
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where  A  is  the  Helmholtz  free  energy  for  n  moles  of  s  species.  This  is  given  by  [19]: 


A  =  Aideal  +  E0(V)  +  nRT/«f(V,T) 


so 


(41) 


_  ideal  +  +  _L  {nRTM(V, T)} 

3nj  dn. 


(42) 


The  ideal  part  of  the  chemical  potential  can  be  taken  as  [1]: 

Hi  (ideal)  =  =  (F°  -  Hg)j  +  (H°)j  +  RTtel 

anj 


f  XjP^ 

vPoy 


(43) 


The  second  term  in  the  expression  for  the  chemical  potential  can  be  written  as: 


dE0  _  dE0  de0  dE0  dV 
dn*  de0  dn,  dn; 


(44) 


where 


27e0 

5 


i  r  i 


f  *  \ 

— 

V 

3 

l-l 

V 

M 

exp 

lVJ 

3v 


( 

V 


(45) 


(46) 


de0 

2^  n  jCfj 

.  j 

dnj 

n 

_ 

n 

dV* 

21  njv-j 

....  J 

*  i 
>1 

i 

dnj 

n 

n 

(47) 


(48) 


The  final  term  in  the  expression  for  the  chemical  potential  can  be  written  as: 


—  (nRT/nf)  =nRT-~  +  KTlnf 


where 

df  =  df  de0  |  df  dV*  |  df  dn 
dr\[  de0  dn^  3v*  dnj  dn  dn^ 


(49) 


(50) 
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Let 


Oi  = 


6eCl 

7.5nRT 


and  a2  = 


54e 


57tnRT 


(51) 


then 


df_ 

9en 


=  ai 


f  M 
v 


vv/3y 


^lna^o) 


2  ( _  *\2 


+  a2 


vv/3y 


6(lnaie0) 


,5  y  *  A3 


+  a3 


vW3y 


9(lna!e0)8 


+3o2V°2e0 


f..*V3 

-2 


n 


vvy 


Vv2 


LV 


exp 


n 


y  ^ 


-/ 


v  v  y 


7J 


(52) 


Now, 


3f  _  9f  9a  j  9f  9a2 

9n  9aj  9n  9a2  9n 


(53) 


where 


9a  i  ax 

9n  n 

9 a2  _  a2 
9n  n 


(54) 


(55) 


9f 

y  *  a 

V 

3(lnaieo)2  _ 

(  *  \ 
V 

60naieo)5  ,  _ 

y  *  \ 

V 

6  9(lnaie0)8 

9ai  1 

v/3 

t  ao 

<*i 

lv/3J 

_  +a3 

<*1 

U3J 

(56) 


and 


9v 


/3v 


2 

df  a^lna^)3  2a2v*(lnaxe0)6  3a3|v  )  0noleo) 

(/3y)2  (/3v)3 


+3(a2e0)2 


-t  y  *  v 

V 

j 

V 

-2 

V 

V 

^  y 

l 2 

/  *\-2 

2/ 

v  v  +  — 

3 

\  / 

3' 

n  (  i\ 
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v  v  y 


4  1 


exp 
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The  internal  energy  can  be  derived  using  the  same  approach  adopted  by  Mader  with 
the  BKW  equation  of  state: 


E(T,  v)  =  nE(T,  °°)  +  nRT2  j  ^  f “  Jdv  +  E0  (v) 


(58) 


where 


E(T,~)=Ixi(E^)j 

The  G  factor  can  be  written  as: 

G(T,v)=  1  +  3 
v 


(59) 


(60) 


where 


u  =  c^x3  +  2a2x6  +  3a3x9  -  a6<*-~ 

T2 


(61) 


q  f.  o  CU 

v  =  1  +  oqx  +  a2x  +  a3x*  +  — -■ 

T2 


(62) 


ai  =  ai 


f  *  \ 

V 


,v/3y 


(63) 


a2  =a2 


(  *  \ 

v 


\yV  ) 


(64) 


a3  =  a3 


(  M 
'  v  1 


Vv/'  j 


(65) 


a4  = 
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3 

-2 

(  *\ 

V 

V 

V 

V 

V  ) 

V  y 

x3 


(66) 


a5  =  2x2a4.N/a7 


(67) 
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a6  =^2Vo4 


a?  -  x3 


117 

4 


(68) 


(69) 


6e0eCl 
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54e0e 
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r 
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x  =  ln\ 
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The  integral  in  equation  58  is  evaluated  numerically. 


(73) 


(74) 


(75) 


(76) 


5.3  Code  Validation 

The  equations  outlined  in  the  previous  section  were  programmed  into  SDA*FOR  and 
the  code  was  then  used  to  calculate  the  CJ  state  of  11  representative  CHNO  explosives 
using  the  JCZ3  EOS.  Table  7  shows  the  computed  CJ  pressure,  velocity  and 
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temperature  for  these  explosives,  and  compares  the  results  with  values  obtained  by 
other  authors  using  the  same  equation  of  state  and  different  codes  (primarily  the 

TIGER  code). 

Table  7  shows  that  there  is  considerable  disagreement  between  the  results  obtained  by 
different  authors  using  the  same  equation  of  state.  For  RDX  for  example  t  ree 
calculated  values  of  the  CJ  pressure  vary  from  0.308  Mbar  to  0.322  Mbar,  a  rangeo  °, 
while  for  TNT  four  different  authors  obtain  values  varying  from  0.177  Mbar  to  0.190 
Mbar  for  the  CJ  pressure,  which  is  a  range  of  7%.  Our  computed  CJ  pressure  for 
of  0.311  Mbar  is  within  the  range  obtained  by  other  authors,  and  our  vaue  o  • 
Mbar  for  TNT  is  less  than  0.5%  outside  the  range  quoted  by  others.  For  PETNand  NM 
our  calculated  CJ  pressures  are  also  within  the  range  of  values  calculated  by  other 
authors,  while  for  DATB  and  NQ  our  results  differ  by  no  more !  *an  ^/ofrcim _other 
literature  values.  Our  least  accurate  results  appear  to  be  obtained  for  TATB,  NG  and 
TETRYL,  where  our  computed  CJ  pressures  differ  by  5%  to  7%  from  values  quoted  in 
the  literature.  Even  these  differences  are  still  within  the  limits  of  the  disagreements 
found  by  other  authors  however.  Variations  in  computed  CJ  velocities  are  much  less, 
and  our  values  are  generally  within  2%  of  those  obtained  by  other  authors.  There  is 
much  less  data  on  CJ  temperature,  but  our  computed  values  are  also  in  good 
agreement  with  the  available  data. 

The  values  presented  in  Table  7  suggest  that  the  JCZ3  EOS  has  been  implemented 
correctly  into  our  code  and  we  are  now  able  to  consider  its  application  to  the  CJ  state  of 
PBXN-111.  Before  doing  so  however  it  is  instructive  to  consider  some  of  the  possible 
sources  of  the  disagreements  shown  in  Table  7. 

One  potential  source  of  error  in  the  calculations  is  the  use  of  a  polynomial  fit  to  the 
thermodynamic  functions  of  the  gases  and  solids  (that  is,  F°  and  Hg  in  equation 
(43)).  Historically  this  fit  was  used  to  save  space  and  time  in  the  calculation  but  wit 
the  development  of  computer  hardware  there  is  no  restriction  on  calculating  e 
thermodynamic  functions  exactly.  The  exact  thermodynamic  functions  for  gases 
specified  in  Mader  [1]  (appendix  F)  have  been  incorporated  into  the  program.  It  is  not 
expected  that  the  polynomial  fit  to  the  solid  thermodynamic  functions  will  have  much 
effect  on  the  calculations. 

The  results  for  several  explosives  are  shown  in  Table  8,  where  F  refers  to 
]f°-H°^IT  and  H  refers  to  [h°  -H°).  For  most  of  the  explosives  the  CJ 

pressures  generally  show  only  a  small  deviation  between  the  exact  thermodynamic 
calculation  and  the  approximate  polynomial  fit,  although  for  one  of  the  compositions 
(PETN)  the  difference  is  as  high  as  7%,  and  others  (PETN,  NG,  TET  ) 
differences  of  the  order  of  3%  to  4%.  The  CJ  velocities  show  much  less  variation  than 
the  CJ  pressures,  most  differences  are  on  the  order  of  1%  to  2%,  with  the  larges 
difference  of  2.5%  occurring  for  HNS.  The  CJ  temperatures  show  greater 
discrepancies,  with  PETN  and  NG  both  showing  differences  of  the  order  of  4/o.  Given 
that  differences  of  this  order  can  arise  due  to  differences  in  the  way  the  ideal  quantities 
are  parametrised,  the  differences  observed  in  Table  7  are  less  surprising. 

Another  aspect  of  the  code  which  may  have  had  a  significant  effect  on  the  final  result 
was  the  method  used  to  compute  the  CJ  point.  With  this  in  mind,  we  included  an 
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alternative  algorithm  based  on  the  method  described  by  Mader  (see  [1],  p  442).  In  this 
approach  the  detonation  velocity  is  calculated  at  three  different  pressure  values  and  a 
parabola  is  fitted  to  the  data.  The  minimum  of  the  parabola  then  defines  the  CJ  point. 
Values  calculated  using  this  method  however  showed  no  discernible  difference  from 
the  results  calculated  using  the  previous  method  of  calculating  the  C-J  point,  that  is, 
equation  (23). 

Table  7:  Computed  CJ  Pressure  (Mbar),  Velocity  (m/s)  and  Temperature  (K)for  the  JCZ3  EOS 
from  the  AMRL  code  compared  with  calculated  values  obtained  by  other  authors. 


Explosive 

Density 

AMRLJCZ3 

r  OTHERS  JCZ3 

REFS 

(g/ cm3) 

P 

D 

T 

P 

D 

T 

RDX 

1.8 

0.311 

8741 

4145 

0.322 

8806 

4012 

27 

0.343 

8813 

28 

0.308 

8670 

25 

PETN 

1.77 

0.281 

8220 

4615 

0.288 

8210 

4237 

23 

0.280 

8200 

25 

HMX 

1.89 

0.349 

9133 

3976 

0.356 

9110 

3726 

23 

0.352 

9160 

25 

TNT 

1.64 

0.176 

6786 

3799 

0.181 

6790 

3501 

23 

0.190 

6912 

3647 

24 

0.177 

6911 

3692 

26 

0.188 

6910 

25 

NM 

1.13 

0.121 

6182 

3755 

0.119 

6110 

3467 

23 

0.121 

6245 

3515 

24 

0.117 

6120 

25 

HNS 

1.69 

0.198 

7050 

4300 

0.205 

7140 

25 

TATB 

1.85 

0.254 

7955 

3033 

0.267 

8053 

2957 

24 

DATB 

1.79 

0.240 

7748 

3653 

0.238 

7624 

3269 

24 

0.237 

7700 

25 

TETRYL 

1.70 

0.223 

7476 

4242 

0.240 

7607 

4065 

24 

NQ 

1.55 

0.187 

7384 

2581 

0.191 

7453 

2474 

24 

NG 

1.60 

0.230 

7715 

4927 

0.217 

7535 

4445 

24 

25 
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Table  8:  Computed  CJ  Pressure,  Velocity  and  Temperature  using  Mader  s 
parametrisation  to  the  thermodynamic  functions  and  the  exact  thermodynamic 
functions 


CJ  Pressure  (Mbar) 


Explosive 

Density 

Exact  F 

Exact  F 

Approx.  F 

(g/ cm3) 

Exact  H 

Approx.  H 

Approx.  H 

RDX 

1.8 

0.310 

0.312 

0.311 

1.0 

0.099 

0.103 

0.095 

PETN 

1.77 

0.270 

0.276 

0.281 

HMX 

1.89 

0.350 

0.351 

0.349 

TNT 

1.64 

0.173 

0.173 

0.176 

NM 

1.13 

0.117 

0.117 

0.121 

HNS 

1.5 

0.145 

0.148 

0.156 

TATB 

1.85 

0.257 

0.256 

0.254 

DATB 

1.79 

0.241 

0.241 

0.240 

TETRYL 

1.70 

0.217 

0.219 

0.223 

NO 

1.55 

0.189 

0.189 

0.187 

NG 

1.60 

0.222 

0.229 

0.230 

CJ  Velocity  (m/s) 


Explosive 

Density 

Exact  F 

Exact  F 

Approx.  F 

(g/cm3) 

Exact  H 

Approx.  H 

Approx.  H 

RDX 

1.8 

8709 

8763 

8741 

1.0 

5903 

6017 

6010 

PETN 

1.77 

8106 

8190 

8220 

HMX 

1.89 

9118 

9159 

9133 

TNT 

1.64 

6745 

6775 

6786 

NM 

1.13 

6044 

6084 

6182 

HNS 

- T5 - 1 

6294 

6363 

6454 

TATB 

1.85 

7983 

7982 

7955 

DATB 

1.79 

7744 

7765 

7748 

TETRYL 

1.70 

7403 

7461 

7476 

NO 

1.55 

7422 

7420 

7384 

NG 

1.60 

7642 

7731 

7715 
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CJ  Temperature  (K) 


Explosive 

Density 

Exact  F 

Exact  F 

Approx.  F 

(g/ cm3) 

Exact  H 

Approx.  H 

Approx.  H 

RDX 

1.8 

4046 

4167 

4145 

1.0 

4624 

4836 

4680 

PETN 

1.77 

4432 

4620 

4615 

HMX 

1.89 

3896 

3990 

3976 

TNT 

1.64 

3779 

3834 

3799 

NM 

1.13 

3739 

3788 

3755 

HNS 

1.5 

4313 

4452 

4388 

TATB 

1.85 

3045 

3041 

3033 

DATB 

1.79 

3628 

3672 

3653 

TETRYL 

1.70 

4153 

4277 

4242 

NQ 

1.55 

2591 

2589 

2581 

NG 

1.60 

4743 

4962 

4927 

5.4  Application  to  PBXN-111 

We  now  consider  the  CJ  state  of  PBXN-111  using  the  JCZ3  equation  of  state.  The 
intermolecular  potential  parameters  (ie  n  and  eO  used  for  the  gas  phase  species  are 
listed  in  Table  9.  Values  for  these  parameters  are  the  same  as  those  used  by 
Co wperth waite  and  Zwisler  [19]  with  the  exception  of  HC1  and  CI2.  For  these  two 
molecules  values  of  r,  and  Ei  are  not  available,  so  the  values  shown  in  Table  9  were 
estimated.  The  estimation  was  based  on  the  difference  between  the  intermolecular 
parameters  for  the  exponential-6  potential  and  a  Lennard-Jones  potential.  Values  of  n 
and  £i  for  the  Lennard-Jones  potential  are  available  for  a  large  number  of  species 
including  HC1  and  CI2  [29].  The  average  difference  in  n  and  ei  between  the 
exponential-6  and  Lennard-Jones  potentials  for  a  number  of  molecules  (H2O,  CO2,  N2, 
CO,  H2,  NO,  O2  and  CH4)  was  used  to  adjust  the  Lennard-Jones  parameters  of  HC1 
and  CI2  to  obtain  estimates  of  the  exponential-6  parameters  for  these  molecules.  It 
should  be  noted  that  the  NH3  parameters  shown  in  Table  9  are  assumed  to  be  the 
same  as  those  for  the  isoelectronic  molecule,  H2O.  This  assumption  was  used  by 
Co  wperth  waite  and  Zwisler  [19],  and  it  could  influence  the  amount  of  NH3  in  the 
equilibrium  composition.  The  Cowan-Fickett  equation  of  state  was  used  for  the  solid 
products. 
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Table  9:  Intermolecular  potential  parameters. 


Molecule 

n  (10-10m) 

Ei(K) 

h2o 

3.350 

138.0 

n2 

4.050 

120.0 

co2 

4.200 

200.0 

CO 

4.050 

120.0 

o2 

3.730 

132.0 

h2 

3.340 

37.0 

NHs 

3.350 

138.0 

NO 

3.970 

105.0 

ch4 

4.290 

154.0 

HC1 

3.385 

370.6 

Cl2 

4.495 

367.6 

The  results  of  the  performance  calculations  are  summarised  in  Tables  10  and  11,  along 
with  the  results  obtained  using  the  BKW  equation  of  state,  and  the  ICI IDEX  code  [14], 
The  IDEX  code  uses  the  THEOSTAR  EOS  [30],  which  is  also  based  on  an 
intermolecular  potential.  Table  10  shows  that  both  equations  of  state  based  on 
intermolecular  potentials  provide  a  CJ  state  in  better  agreement  with  the  experimental 
observations,  while  the  parametrised  BKW  EOS  overestimates  both  the  CJ  velocity  and 
pressure  significantly.  There  is  good  agreement  between  the  detonation  velocity 
calculated  using  the  AMRL  implementation  of  the  JCZ3  EOS  (6606  m/ s)  and  the  IDEX 
THEOSTAR  EOS  (6661  m/s),  the  difference  between  them  being  less  than  1%,  while 
the  detonation  pressures  agree  to  within  approximately  10%  (0.201  Mbar  for  JCZ3, 
0.222  Mbar  for  THEOSTAR).  Clearly,  there  are  significant  differences  in  the 
composition  (and  temperature)  calculated  by  each  EOS,  but  there  is  no  experimental 
data  available  to  assess  the  accuracy  of  the  different  compositions. 


Table  10:  CJ  parameters  for  PBXN -111  calculated  with  different  equations  of  state. 


Detonation  Parameter 

BKW 

IDEX 

AMRL 

Pressure  (Mbar) 

0.294 

0.222 

0.201 

Velocity  (m/s) 

8153 

6661 

6606 

Temperature  (K) 

5229 

5297 

5864 

While  the  agreement  between  the  detonation  velocity  and  pressure  calculated  from  the 
two  different  equations  of  state  based  on  intermolecular  potentials  is  reasonable,  there 
is  still  a  significant  difference  from  the  experimentally  determined  values.  Forbes  et  al. 
have  reported  the  detonation  pressure  of  PBXN-111  as  0.120  Mbar,  with  an  infinite 
diameter  detonation  velocity  of  6195  m/s  [5],  while  Bocksteiner  et  al.  have  reported 
the  infinite  diameter  detonation  velocity  of  (Australian)  PBXN-111  to  be  5650  m/s. 
Comparison  between  calculated  and  measured  detonation  parameters  for  highly  non¬ 
ideal  explosives  requires  considerable  care  however.  The  calculated  detonation 
velocities  should  be  compared  with  experimental  values  determined  on  cylindrical 
charges  having  infinite  diameter.  In  practice,  the  experiments  have  been  carried  out  on 
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charges  having  diameters  no  greater  than  twice  their  respective  failure  diameters.  The 
infinite  diameter  results  have  then  been  found  by  extrapolation  using  either  a  liner  [5] 
or  elliptic  [8]  fitting  procedure.  This  is  a  relatively  small  diameter  range  for  non-ideal 
explosives  such  as  PBXN-111  however,  and  experience  with  commercial  composite 
explosives  has  shown  that  diameters  up  to  twenty  times  the  failure  diameter  need  to 
be  used  before  reliable  estimates  of  the  infinite  diameter  detonation  velocity  can  be 
made  [31]. 


Table  11:  CJ  composition  (moles  per  100  g  of  explosive)  for  PBXN-111  calculated  with 
different  equations  of  state. 


Product 

BKW 

IDEX 

TCZ3 

H,0 

0.6121 

0.2681 

0.3620 

H, 

0.6787 

0.2775 

0.0432 

cm 

0.0877 

0.1556 

0.0298 

co7 

0.0019 

0.0067 

0.0883 

CO 

0.0455 

0.3436 

0.1227 

n7 

0.4239 

0.2197 

0.1055 

NH, 

0.0633 

0.4661 

0.7000 

NO 

- 

0.0009 

- 

o7 

0.0000 

0.0000 

0.0001 

HC1 

0.2839 

0.2591 

0.3650 

Cl, 

0.0411 

0.0534 

0.0005 

C(graphite) 

0.9478 

0.0000 

0.8393 

C(diamond) 

- 

0.6338 

- 

ai7o^ 

0.4635 

0.4632 

0.4635 

The  detonation  pressure  reported  by  Forbes  is  based  on  measurements  of  the 
underwater  shock  velocity  obtained  from  a  standard  aquarium  test  conducted  on 
PBXN-111  [32].  The  pressure  is  then  calculated  from  the  velocity  measurement  using  a 
procedure  which  was  derived  for  explosives  having  ideal  detonation  behaviour. 
Forbes  notes  that  this  procedure  is  questionable,  and  may  result  in  errors  of  up  to  30% 
in  the  reported  detonation  pressure.  Kennedy  and  Jones  [14],  and  more  recently 
Kennedy  [33],  have  discussed  the  behaviour  of  non-ideal  explosives  such  as  PBXN-111 
in  more  detail,  and  have  shown  that  the  concept  of  a  unique  "detonation  pressure"  (or 
C-J  pressure)  has  little  significance  for  such  explosives. 

The  detonation  parameters  for  PBXN-111  recorded  in  Table  10  were  calculated  using 
the  method  of  steepest  descent  algorithm.  Use  of  the  probabalistic  algorithm  would 
have  required  considerably  more  computational  time,  and  in  view  of  the  results 
obtained  in  Section  4,  and  the  accuracy  of  some  of  the  potential  parameters  used  in  the 
calculation,  this  was  considered  to  be  unnecessary. 
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6.  Discussion  and  Conclusion 

This  report  has  described  the  development  and  application  of  a  chemical  equilibrium 
computer  code  to  calculate  the  detonation  parameters  of  condensed  phase  explosives. 
The  code  is  based  on  the  Chapman-Jouguet  theory  of  detonation,  incorporates  two 
different  equations  of  state  to  describe  the  detonation  products,  and  offers  a  choice  of 
two  different  algorithms  to  minimize  the  free  energy  of  the  product  composition. 

Using  the  BKW  EOS  and  either  the  method  of  steepest  descent  or  the  probabalistic 
algorithm  to  minimise  the  free  energy  the  code  accurately  reproduces  the  detonation 
parameters  of  ideal  CHNO  explosives.  The  method  of  steepest  descent  is  inherently 
faster,  but  is  not  guaranteed  to  find  the  global  minimum.  The  probabalistic  algorithm 
is  considerably  slower  but  has  the  advantage  of  not  requiring  an  initial  estimate  of  the 
equilibrium  composition,  and  is  guaranteed  to  find  the  global  minimum.  Either 
minimization  method  can  be  used,  depending  on  the  particular  application,  and  the 
amount  of  computer  time  available. 

As  well  as  BKW  the  code  also  employs  the  JCZ3  equation  of  state.  BKW  contains  four 
adjustable  parameters,  and  these  have  been  optimized  to  provide  good  agreement 
with  detonation  parameters  for  ideal  CHNO  explosives.  JCZ3  is  based  on  an 
intermolecular  potential  and  contains  no  adjustable  parameters,  other  than  the  well 
depth  and  equilibrium  bond  distance  for  each  interaction.  JCZ3  was  found  to  be 
considerably  better  than  BKW  in  predicting  the  detonation  parameters  of  the 
underwater  explosive  PBXN-111,  and  we  expect  that  in  general  JCZ3  will  be  the  more 
appropriate  equation  of  state  to  use  for  highly  non-ideal  explosives. 

The  results  obtained  here  show  that  considerable  care  should  be  exercised  when 
attempting  to  use  the  code  to  calculate  detonation  parameters  for  a  particular 
explosive.  The  code  is  based  on  the  Chapman-Jouguet  theory  of  detonations,  and  as 
such  is  only  applicable  to  explosives  which  behave  ideally.  When  used  to  model 
explosives  which  behave  non-ideally  the  predicted  detonation  velocity  and  pressure 
require  careful  interpretation. 
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9.  Appendix:  Running  SDA.FOR  and  PEA.FOR 


This  appendix  describes  the  actual  usage  of  the  two  new  programs  for  computing  CJ 
properties  of  explosives,  SDA.FOR  and  PEA.FOR.  Firstly  the  input  file  for  these  codes 
will  be  described  followed  by  a  description  of  the  output  that  each  of  the  codes 
produces.  The  final  section  will  reveal  a  number  of  useful  hints  for  achieving 
successful  results  from  these  codes. 


INPUT  DESCRIPTION 

Since  the  input  files  for  SDA.FOR  and  PEA.FOR  are  virtually  identical  we  will  only 
describe  the  file  once  with  a  clear  indication  of  the  minor  differences  where  applicable. 
The  input  file  is  mostly  free  format  and  Table  12  describes  each  line  of  input  in  detail. 

There  are  two  differences  in  this  input  for  the  SDA  and  PEA  programs.  In  card  8  and 
card  15  an  initial  estimate  of  the  number  of  moles  of  each  product  is  entered.  Both 
SDA  and  PEA  read  in  this  input  but  only  the  SDA  program  actually  uses  it.  The  PEA 
program  ignores  initial  guesses.  The  second  difference  occurs  in  card  17.  Only  the 
PEA  program  requires  this  input  card,  the  SDA  program  does  not  read  in  the  nlim 
value. 


OUTPUT  DESCRIPTION 

Table  13  shows  a  sample  of  the  output  from  a  run  of  the  SDA.FOR  program.  The 
output  of  the  PEA.FOR  program  is  identical.  The  output  below  represents  the  last 
three  cycles  in  the  calculation  of  the  detonation  parameters  of  pure  RDX  at  a  loading 
density  of  1.8  g/cm3.  Each  cycle  represents  a  further  refinement  of  P  and  T  as 
described  earlier  in  the  four  step  iterative  process  for  determining  the  CJ  point  (see  just 
below  equation  28).  The  output  at  each  cycle  prints  the  current  value  of  the 
temperature  and  pressure  followed  by  a  list  of  each  product,  the  current  number  of 
moles  of  that  product  and  the  amount  each  product  was  changed  by  during  the  last 
step  of  the  equilibrium  composition  determination.  The  last  line  of  output  of  each 
cycle  reprints  the  pressure  and  a  value  called  fvp.  The  fvp  value  is  a  measure  of  how 
close  equation  28  actually  is  to  being  zero.  When  fvp  falls  below  a  specified  tolerance 
then  equation  28  is  deemed  to  be  satisfied  and  the  program  will  terminate.  Prior  to 
terminating,  the  program  calculates  the  detonation  velocity  and  finally  prints  out  the 
detonation  temperature,  pressure  and  velocity  as  shown  at  the  very  bottom  of  the 
sample  output. 


RUNNING  HINTS 

There  are  a  number  of  points  that  should  be  kept  in  mind  when  using  the  SDA.FOR 
and  the  PEA.FOR  programs.  Experience  has  shown  that  under  certain  circumstances 
the  program(s)  will  fail  to  give  a  successful  result  unless  these  points  are  carefully 
considered  when  preparing  the  input  file.  In  most  cases  when  the  program  fails  to 
produce  a  result,  the  output  it  does  produce  usually  indicates  the  nature  of  the 
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problem  and  the  following  hints  describe  the  approach  that  has  proved  successful  in 
overcoming  the  original  problem. 

1)  The  fvp  value  is  diverging  away  from  zero. 

Under  certain  circumstances  the  program  will  adjust  P  in  such  a  way  that  the  fvp 
value  does  not  tend  to  zero  as  it  should.  This  usually  happens  if  the  initial  guess  of 
the  detonation  pressure  is  very  poor  and  the  problem  seems  to  be  worse  for  explosives 
with  a  low  loading  density.  The  remedy  is  simple,  restart  the  calculation  with  a 
different  estimate  of  the  detonation  pressure. 

2)  Solid  products  with  a  large  free  energy. 

If  any  of  the  solid  products  has  a  large  free  energy  then  the  equilibrium  composition  of 
that  product  will  tend  towards  zero.  If  the  equilibrium  composition  of  a  solid  product 
becomes  zero  the  calculation  will  fail.  This  will  be  evident  as  the  program  will  print 
the  word  "solid"  followed  by  a  number  which  indicates  which  solid  is  causing  the 
trouble.  This  is  followed  by  a  message  which  reads  "rlam  =  0".  If  this  occurs  the  only 
solution  is  to  remove  that  solid  from  the  product  list  and  restart  the  calculation.  In 
general  it  may  be  prudent  to  start  a  calculation  with  no  solid  products  and  then  add 
each  solid  product  one  at  a  time  to  see  if  any  of  the  solid  products  is  going  to  cause  a 
problem. 

3)  Initial  guess  of  equilibrium  composition. 

This  is  only  relevant  for  the  SDA.FOR  program  since  the  PEA. FOR  program  does  not 
require  an  initial  estimate  of  the  equilibrium  composition.  In  general  it  is  best  to 
assume  that,  for  CHNO  explosives,  the  major  detonation  products  will  be  H20,  C02, 
CO  and  N2.  The  initial  guess  of  the  equilibrium  composition  should  then  comprise 
appropriate  amounts  of  these  products.  Thus  all  of  the  H  in  the  explosive  formulation 
will  become  H20,  all  of  the  N  will  become  N2  and  the  amounts  of  CO  and  C02  will  be 
determined  by  the  oxygen  balance  of  the  explosive  formulation.  If  the  explosive  is 
particularly  oxygen  deficient  then  more  CO  will  be  required  to  maintain  mass  balance. 
Other  gaseous  products  can  still  be  included  in  the  list  of  possible  products  and  can  be 
given  a  very  small  initial  value  for  y(i).  Solid  products  can  also  be  included  with  a 
small  initial  value  of  y(i).  It  is  not  necessary  to  have  an  initial  guess  that  satisfies  mass 
balance  although  it  is  suggested  that  the  initial  guess  should  not  be  too  far  away  from 
mass  balance  to  avoid  potential  problems  in  the  determination  of  the  correct 
equilibrium  composition. 
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Table  12: 


Line 

Variable 

Name(s) 

Type 

Format 

Description 

1 

nelem 

integer 

free 

number  of  unique  elements  in  explosive 
formulation 

2 

react(i) 

real 

free 

for  each  element  in  the  explosive  formulation  this 
array  holds  the  amount  of  that  element  present 

3 

hfe,rmme, 

dens 

real 

free 

hfe  is  the  explosive  formulation  heat  of  formation  at 
OK  in  cal /formula  weight,  rmme  is  the  formula 
weight  of  the  explosive  and  dens  is  the  density  of 
the  explosive. 

4 

ngprod 

integer 

free 

number  of  gaseous  detonation  products 

5 

gases(i) 

string 

a5 

label  description  for  gaseous  product 

6 

prodg(i,j),  j=l  to 
nelem 

real 

free 

specifies  amount  of  each  element  in  this  gaseous 
product 

7 

ag(i),bg(i), 
cg(i),dg(i), 
eg(i),ricg(i),  hf(i) 

real 

free 

ag,bg,cg,dg,eg  and  ricg  are  the  constants  used  in 
equation  8  and  hf  is  the  heat  of  formation  of  the 
gaseous  product  at  0  K. 

8 

y(i),rk(i) 

real; 

free 

y  is  the  initial  guess  of  the  number  of  moles  of  the 
gaseous  product  and  rk  is  the  co-volume  of  the 
gaseous  product  (see  equations  2  and  3) 

Cards  5  to  8  are  repeated  for  each  gas  species  (ie  ngi 

a  rod  times) 

9 

nsprod 

integer 

free 

number  of  solid  detonation  products 

10 

solids(i) 

__strin8 

a5 

label  description  for  solid  product 

11 

prods(i,j),  ]-l  to 
nelem 

real 

free 

specifies  amount  of  each  element  in  this  solid 
product 

12 

as(i).bs(i), 

cs(i),ds(i), 

es(i) 

real 

free 

these  are  the  Cowan  solid  equation  of  state 
parameters  (see  for  example  equation  10) 

13 

als(i),a2s(i),cls(i) 

,c2s(i),c3s(i) 

real 

free 

these  are  the  remaining  Cowan  solid  equation  of 
state  parameters 

14 

ase(i),bse(i),cse(i) 

,dse(i),ese(i),rics( 

i),hfs(i) 

real 

free 

ase,bse,cse,dse,ese  and  rics  are  the  solid  product 
equivalent  of  the  constants  used  in  equation  8  and 
hfs  is  the  heat  of  formation  of  the  solid  product  at  0 
K. 

15 

ys(i),rmm(i),vo(i) 

real 

free 

y  is  the  initial  guess  of  the  number  of  moles  of  the 
solid  product,  rmm  is  the  molecular  mass  of  the 
solid  product  and  vo  is  the  initial  volume  (in 
cm3/g). 

Cards  10  to  15  are  repeated  for  each  solid  species  (ie  nsprod  times) 

16 

press, temp 

real 

free 

starting  values  for  the  pressure  and  temperature 

17 

nlim 

real 

free 

number  of  guesses  used  in  the  probabilistic 
algorithm,  this  input  is  only  required  for  the 
PEA.FOR  program. 
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Table  13: 


Tcj= 

2587.792968750000 

Pcj= 

0.3451171875000000 

H20 

2.998805060393444 

2.0324615110367539E-08 

C02 

1.489432550814017 

1.6185742568453065E-07 

N2 

2.999984618027165 

2.9389973510429712E-10 

H2 

1 .1 487936880493267E-03 

-1 .9442915771799424E-08 

CO 

2.2324121316814211E-02 

-3.4390013300231604E-07 

NH3 

3.0763945670919224E-05 

-5.8779959834121635E-10 

02 

2.8583308535342730E-06 

-6.9666778548661276E-1 1 

SOLC 

1.488243327869169 

1.8204270731819738E-07 

P 

0.3451171875000000 

fvp 

2.9294807983619442E-04 

Tcj= 

2587.792968750000 

Pcj= 

0.3450976562500000 

H20 

2.998804735844558 

2.0336275741161813E-08 

C02 

1.489430558980391 

1. 61931 95202102162E-07 

N2 

2.999984612673078 

2.9410318902467480E-10 

H2 

1 .1491021 746754503E-03 

-1.9453966329874089E-08 

CO 

2.2328430423178673E-02 

-3.4406089971365274E-07 

NH3 

3.0774653844636130E-05 

-5.8820624758797982E-10 

02 

2.8578857408429585E-06 

-6.96399882034161 70E-1 1 

SOLC 

1.488241010596430 

1.8212894770663901E-07 

P 

0.3450976562500000 

fvp 

2.8134827164343762E-04 

Tcj= 

2587.792968750000 

Pcj- 

0.3450878906250000 

H20 

2.998804573543107 

2.0444491077853400E-08 

C02 

1.489429562965963 

1 .6268356730908540E-07 

N2 

2.999984609995493 

2.9605989049219516E-10 

H2 

1 .1492564433704240E-03 

-1 .955631 1564429096E -08 

CO 

2.2330585198497803E-02 

-3.4567291897485020E-07 

NH3 

3.0780009014973331 E-05 

-5.9211970592613272E-10 

02 

2.8576632343753382E-06 

-6.93534588031 25902E-1 1 

SOLC 

1.488239851835539 

1. 82989351 64640096E-07 

P 

0.3450878906250000 

fvp 

1. 2516557871 852102E-04 

Tcj= 

Dcj= 

2587.792968750000  K 
8711.452096205158  m/s 

Pcj= 

0.3450781250000000  Mbar 
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Amount  of  AP  reacting 
(g  per  lOOg  of  explosive) 
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Figure  2:  Effect  of  extent  of  reaction  of  ammonium  perchlorate  (AP)  on  the  calculated 
performance  of  PBXN-1 1 1 . 
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20.  ABSTRACT 


We  present  a  detailed  description  of  the  development,  implementation,  and  application  of  a  computer  program  to  calculate  the 
detonation  parameters  of  condensed  phase  explosives.  The  code  is  based  on  Mader's  BKW  chemical  equilibrium  code,  but 
contains  important  new  features.  A  new  algorithm  to  calculate  the  minimum  in  the  free  energy  of  the  product  composition  has 
been  included.  This  is  a  probabilistic  algorithm,  based  on  the  method  of  Benke  and  Skinner,  and  its  inclusion  ensures  that  the  true 
global  minimum  in  the  free  energy  will  always  be  found.  As  well  as  the  BKW  equation  of  state  to  describe  the  detonation 
products,  the  new  code  also  includes  the  JCZ3  equation  of  state.  This  is  an  intermolecular  equation  of  state  containing  no 
adjustable  parameters,  and  hence  should  be  applicable  to  a  wider  range  of  explosives  than  could  be  considered  using  the  BKW 
code.  We  have  validated  the  code  on  a  wide  range  of  military  explosives,  using  both  the  new  probabilistic  minimisation  algorithm 
as  well  as  the  original  method  of  steepest  descent,  for  both  the  BKW  and  JCZ3  equations  of  state.  We  also  present  a  detailed 
description  of  the  application  of  the  code  to  the  non-ideal  underwater  explosive  PBXN-111,  and  show  that  the  performance  of  the 
explosive  is  best  described  using  the  JCZ3  equation  of  state. 


